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1. Introduction 

Metamaterials have large potential in electromagnetics and optics due to their possibility of 
tailoring the permittivity and permeability. This enables construction of, for example, media 
with negative refractive index n [1], perfect lenses [2], invisibility cloaks [3,4], and other 
exciting components transforming the electromagnetic field [5]. Unfortunately, the perfor- 
mance of such devices is strongly limited by losses. Although causality and passivity do 



not prohibit negative index materials with arbitrary low losses [6], in practice it is difficult 
to fabricate materials with high figure of merit FOM = — Ren/Imn, especially at optical 
frequencies [7-10]. For perfect lenses, losses limit the amplification of evanescent waves as- 
sociated with large spatial frequencies, and therefore the resolution [11,12]. It has therefore 
been suggested to introduce gain into the metamaterials [7,13-21]. This could be a promising 
approach provided the intrinsic losses can be made relatively small so that compensation by 
a realistic amount of gain is possible. 

Both permittivity e and permeability /x may involve losses; thus in general, gain may be 
needed to reduce both Ime and Im/i. For a perfect lens it is generally not sufficient e.g. to 
reduce Ime below zero such that the refractive index n = y/ejl becomes real. However, as 
long as the object to be imaged is one-dimensional, only one polarization (TE or TM) of the 
electromagnetic field is required. Then, provided the lens is sufficiently thin, only one of the 
parameters e and fi is relevant for the transmission of evanescent waves [2,12]. Choosing TM 
polarization, only e matters, enabling gain compensation with dielectric, active media. 

Introducing the necessary active material into a metamaterial leads to a change not only 
in the imaginary part but also the real part of the permittivity, and should therefore be kept 
in mind while designing the metamaterial structure. Other critical considerations include 
matching of the negative refractive index frequency band to that of the gain lineshape func- 
tion, the level of loss possible to overcome in the absence of saturation, and the saturation 
constant of the active medium. 

There have been several attempts to create a perfect lens in the near IR-spectrum the 
last years [7-9,22]. The FOM currently reported is of the order of 3 for the frequency where 
Ren ^ —1 [8,9]. With these values, traditional optical amplifiers such as Erbium-doped silica 
or gas laser amplifiers will not be able to reduce the intrinsic losses significantly. Theoretical 
studies have shown that it may be possible to raise the FOM at near IR-frequencies to as 
much as 20, while keeping Ren ^ —1 [23]. It has also been reported that laser dyes, or 
dye-Ag aggregate mixtures, may reach amplifications of up to Imn ^ —0.06 at near IR- 
frequencies [14,24]. Taking into account these reports, this article will not speculate further 
on the choice and design of the metamaterial, but merely assume an appropriate material is 
physically feasible. 

The main purpose of our work is to consider the transmission of evanescent waves in a 
practical, gain-compensated perfect lens. Clearly, gain saturation is highly relevant in this 
context, and we demonstrate how this effect leads to limited amplification of evanescent fields, 
and therefore limited resolution. We calculate the resolution as a function of the saturation 
constant of the active medium, and also the detailed field profile and refiections from the 
lens. It will become clear that gain saturation is a critical effect which may lead to severe 
limitations. 



2. Nonlinear gain saturation and field calculations 

The relative permittivity of the active metamaterial is given by: 

e(a;) = l + XpM + Xa(w), (1) 

where Xp{^) denotes the susceptibility of the passive structure, and Xa{^) the contribution 
from the active part. The contribution to the susceptibility for travelling and evanescent fields 
in active media can be modeled using semiclassical theory [25] . If there is spherical symmetry, 
allowing for coupling to several degenerate states with different values of the quantum number 
m, one can show using the Wigner-Eckart theorem and the general properties of the Clebsch- 
Gordan coefficients that the system can effectively be treated as a two-level system. Assuming 
two-level atoms and using the dipole approximation, we find the following expression for the 
active susceptibility: 



Xa(u;) = ^ ' ' . (2) 

Here Uq is the transition frequency, u the frequency of the incident light, 7 a phenomeno- 
logical decay rate due to spontaneous emission and elastic collisions, E the complex electric 
field, and Ss{u}) the saturation constant of the active medium. The saturation constant de- 
pends on the selected gain material and pumping level. For dye amplifiers a normal value 
is Ss{ujq) ~ 10^ V/m [24,26], corresponding to intensities in the kW/cm^ regime for propa- 
gating waves. The numerator in Eq. (2) describes the susceptibility in the limit |E|/£^s — )■ 0. 
The numerator contains the line shape function and several material parameters; factors 
that are irrelevant for the analysis below are absorbed into the function A{u). For u = Uq 
and |E|/£^s ~^ 0, A^u) is simply —lm.Xa.{^)- Both functions A{uj) and Ss{uj) are real-valued. 
Note that the material parameters describing the active medium are effective parameters, 
depending on the geometry of the metamaterial structure, and are not necessarily equal to 
the bulk parameters of the gain material [18, 19]. 

Throughout this paper, we will consider the frequency u = Uq. Eq. (2) now reduces to: 

-iAjup) 
Xa(wo) = ^^- 3 

Eq. (3) describes how the imaginary part of the total permittivity relates to the pumping and 
local field amplitude. The real part of the total permittivity is independent of the pumping 
and the local field. 

Note that there is a fundamental difference between the nonlinearity due to gain saturation, 
and conventional second- and third-order nonlinearities. First, the nonlinearity due to the 
denominator in Eq. (2) is so large that a Taylor expansion up third order is generally not valid. 



Second, the nonlinearity of a gain medium is expressed in terms of the slowly varying field 
envelope |E|, while the second- and third order nonlinearities usually are expressed in terms 
of the rapidly varying time-dependent electric field. Since the nonlinearity in our case can be 
characterized using |E|, the medium will not generate new frequencies for monochromatic 
input [27]. Nevertheless, if several modes (or frequencies) are present, the modes interact in 
the sense that the complex refractive index seen by one mode is dependent on the presence 
and strength of all modes. 

We consider a perfect lens slab which extends to infinity in the xy-plane, and has thickness 
d in the ^-direction, see Fig. 1. The source is located a distance a (with a < d) from the 
input end of the lens. The incident field from the source will be taken to be a superposition 
of plane TM-waves, with the magnetic field in the y-direction. Provided uod/c <^ 1 and 
\fi\ ~ 1, the spesific value of /i is not critical for the operation of the lens for evanescent 
waves [2,12]. Here c is the vacuum light velocity. The permittivity is given by Eqs. (1), 
(3), and ReXp(wo) = ~2. The remaining losses after gain compensation (in the absence of 
saturation) is described by the parameter: 

Ax = lraXpiuJo)-A{uJo). (4) 

To find the steady state solution to Maxwells equations for our nonlinear medium, an 
iterative approach can be used. In the zeroth iteration, the electric field is simply set to zero 
everywhere. (Alternatively, the initial field could be set to infinity. This does not give any 
significant difference in performance, in terms of the required number of iterations.) In the 
next iteration, Eqs. (1) and (3) are used to find an approximation of the permittivity of 
the lens. Taking the incident magnetic field to be unity (normalized), we can compute the 
magnetic and electric fields everywhere. Now we may repeat the iteration; calculate a new 
approximation of the permittivity from the field, computing the resulting field from this new 
structure, etc. The iteration procedure has an inherent stability, as growing fields leads to 
less gain in the medium, and vice versa. 

Nevertheless, inaccuracies and even divergence may arise if the number of slices is too 
low, so that the field no longer can be treated as constant in each slice. In the case Ax = 
the convergence seems to be particularly sensitive to the number of slices. An alternative to 
increasing the number of slices to a very high number, is to regularize the iterative approach 
as follows: Rather than setting the permittivity to that resulting from the field in the previous 
iteration, it can be set to a weighted mean of the permittivities as resulting from the last 
two iterations. In our computations, the permittivity in iteration i (for i > 2) was set to 
0.5 times the permittivity calculated by the field from iteration i — 1, pluss 0.5 times that 
resulting from iteration i — 2. For i = 1 the permittivity was calculated using the field from 
iteration 0. This resulted in convergence after ~ 20 — 35 iterations. The weight factor 0.5 



is somewhat arbitrary; other choices are possible but may require a larger A^ or number of 
iterations to obtain convergence. 

If the root mean square deviation of three successive iterations are within a specified limit 
(10~^^ for the relative permittivity in our computations), and strictly decreasing, the results 
are deemed converge. Note that when the fields of subsequent iterations coincide, we have 
a valid solution to Maxwell's equations with constitute relation as implied by Eqs. (1) and 
(3). 

In general, the fields in one iteration, and therefore the permittivity in the next iteration, 
will be dependent on x and z. Thus the computation of the fields in the next iteration 
requires the solution to Maxwell's equations in an inhomogeneous structure. Note that in 
each iteration, the structure is linear; the nonlinearity of the structure enters through the 
iteration. For the linear calculation, we employ a transfer matrix technique, considering the 
different plane waves in the structure. The lens is divided into A^ slices in the x|/-plane, as 
seen in Fig. 1. These slices must be sufficiently thin, such that the permittivity inside each 
slice is approximately uniform in the z-direction. For this condition to be valid for the next 
iteration as well, the resulting field from the present iteration must also be approximately 
constant. This means that k^d/N < 1 for the transverse wavenumbers k^ that contribute 
significantly to the fields. 

The electric field can be found using the Ampere-Maxwell's law: 

E(x, z) = — — — V X H(x, z). (5) 

—zue[x, z)eQ 

With periodic boundary conditions in the x-direction, the magnetic field and the permittivity 
can be expanded in discrete Fourier series 

ii{x,z) = H{x,z)y = '^hm{z)exp{ik^mx)y, (6) 

m 

e{x,z) = y^^em{z)exp{ikr,mx), (7) 

m 

for some Fourier coefficients hm{z) and em{z). Here kxm = 'i'Km/L, L is the computational 
domain, and y is the unit vector in the ^/-direction. 

From Maxwell's equations, we find that the magnetic field satisfies 

V^i7 + 6/ifc^i7--|^^ = 0, (8) 

e ox ox 

where k = uq/c. We express 1/e and {l/e)de/dx as Fourier series as follows 

-—- = \^Q^ exp{ikxrnx), (9) 

ex -^^ 



and 

1 de(x) ^-^ ^ , , . 

-r^^- — = > K„exp(zfc^mx). (10) 

eix) ox ^-^ 

The Fourier coefficients Fm are now given as 

m' 

or as a matrix product 

F = iQkxe, (12) 

where F = {-Fm,}™,, kx = diag(/i;2.m), e = {em}m, and Q is a Toephtz matrix with elements 
Inserting the Fourier series into equation (8), we obtain 



dz2 



klm^miz) + k'^^l^em~rn'{z)h^i{z) -^iFm-m'Km'hm' = 0, (13) 



for each m. Let h = {hm}m& and use fc^^ = k"^ — k^.^. We can write Eq. (13) as a matrix 
equation 

^ + [k.^ + V] h(.) = 0, (14) 

where kz = diag(fc2m), and V is the operator defined as 

V = -A;^! + A;VG - F. (15) 

Here G and F are infinite dimensional Toeplitz matrices with elements Gjj- = Ci-j and 

Eq. (14) may be decomposed into a first order system by writing h = h+ + h^, where 
h+ = {h^}m£Z and h~ = {h^}mez- In fact, the equations 

dh+ 

-^ = tKh+ + t{2K)-W{h+ + h-), (16a) 

dh" 

-— = -^k,h"-2(2k,)-iV(h+ + h-), (16b) 

dz 

are seen to be equivalent to (14) after differentiation and summation. This decomposition is 
particularly convenient, since outside the lens V = 0, and Eq. (16) has the simple solution 
h^{z) = const ■ exp{±ikzmz)- In other words, outside the lens, h^ and h^ are the forward 
and backward propagating waves, respectively. 
By adding Eqs (16a) and (16b), we obtain 

^ = 2k,(h+-h-), (17) 



which is needed for V x H. 
By writing 



^-. 



h- 



iK + i{2K)~^v 



-^(2k, 



'V 






(1^ 



Eq. (16) may be brought into matrix form: 

d^ 

d7 



c^. 



(19) 



Since the permittivity is assumed to be independent of z within a shce of the lens, the 
matrix C will be constant for each slice. Thus, Eq. (19) can be integrated to obtain 



^(^fe) = eXjp{{Zh - Za)C}'^{Za) 



(20) 



for Za and Zb inside the same slice. Let Zj be at the left-hand side of slice j, Aj the thickness 
of the slice, and Cj the C- matrix for layer j. From the field at Zj, we find the field at the 
right hand side of the slice as 



*(^j + A,)=exp{A,C,}^(z,). 



(21) 



Note that j = corresponds to the region between the source and the lens, j = 1, ...,N are 
the slices inside the lens, and j = N + 1 represents the region from the lens to the image 

b. Let us define 



plane. The thicknesses are Aq = a, Aj = d/N, for j = 1, ..., A^, and A^v+i 

Mj = expjAjCj}. 



(22) 



Eq. (22) propagates the field from the start of a slice, to the end. Next, connecting the fields 
of adjacent slices with the electromagnetic boundary conditions, we find for the boundary 
between slice j and j + 1 that 



dHj+i{zj+i) _ ej+idHj{zj+i] 



dz 



dz 



(23) 



Inserting the Fourier series from Eqs. (6) and (7), we obtain a convolution on the right hand 
side corresponding to the multiplication of a Toeplitz matrix and the vector containing the 
components dhm{z) /dz. The Toeplitz matrix is defined by the Fourier components of ej+i/e^ 
and will be called P. Then 



dh, 
dz 



dh 



Ui+iJ 



jK'^j+i, 



■i+i 



dz 



.^i+iJ 



(24) 



where P,- corresponds to the transition between layer j and j + 1. 



Eq. (24) together with the fact that h is continuous across the layer boundary, gives us 
the following transfer matrix 






kz rjkz i 
Kz r j Kz ^ 



kz rj kz 

kz r j kz 



(25) 



Let us call the transfer matrix in Eq. (25) Vj. 

By the successive application of (22) and (25), we find: 



T 




M 



N+l 



UiV.M, 



j=N 



I 

R 



(26) 



Here, each column i of I corresponds to an experiment where the incident field amplitude 
is 1 for one of the Fourier components and zero for the others. The z'th column of R is the 
reflection at the source plane of experiment i, and the z'th column of T is the corresponding 
transmission at the image plane. To get the reflection in the case of two or more waves, the 
corresponding columns of R are added, the new transmission are found by adding colums of 
T. Once the total matrix in Eq. (26) has been found, it is straightforward to calculate the 
unknowns T and R, and therefore the field amplitudes in all slices. 

3. Numerical results 

The thickness d of the lens was chosen such that u^d/c = 2tc/10. The resolution clearly 
improves with decreasing distance b from the lens to the image, since then the required 
evanescent fields at the end of the lens are reduced. However there may be practical reasons 
that makes it impossible to reduce the distances a and b below a certain value. In our 
simulations we have taken a = b = d/2. For simplicity we normalized eo = /io = "^o = 1- 
The permeability was set to /i = —1; however, since the lens was relatively thin, the specific 
value of /i did not matter significantly for evanescent TM waves. The number of slices were 
taken to be A^ = 20. The computation domain L was chosen in the range (15, 50) depending 
on the specific problem, and the number of Fourier components (m values) was of the order 
of 100. 

First, a single mode source H = exp(z/cj.x + ikzz)y was considered. The reflection and 
transmission coefficients, and the fields in the lens, were computed using the iterative method 
above. The transmission coefficient is shown in Fig. 2. It is easy to see improvements as a 
result of gain compensation, dependent on the saturation constant £s- (Here it is useful 
to recall that for a normalized magnetic field H = e'XY){ikxX + ikzz)y, the electric field is 
E = {\kz/ujQeeQ\k — [/c^./a;oeeo]z) ex^iiik^x + ikzz). Thus, for our normalization, we see that 
|E| ~ /i;^,., to be compared to Eg.) Nevertheless, for a fixed amplitude of the incident field. 
Fig. 3 indicates that an exponential increase in the saturation constant is needed for a linear 



increase in the resolution. This is an important result, as it shows the difficulty of getting 
large resolution: The required, large evanescent fields associated with large spatial frequencies 
saturate the gain at the output of the lens. 

The refiection coefficient is plotted in Fig. 4. We note that significant refiections arise 
even for the spatial frequencies where the transmission is relatively large. As can be seen in 
Fig. 5, the field distributions of the two evanescent components in the lens increase roughly 
exponentially with +z or —z, respectively. For small spatial frequencies, where the lens 
is essentially perfect, the one increasing in the +z direction dominates. For higher spatial 
frequencies the two components have a similar amplitude, such that the total field and 
therefore the imaginary part of the permittivity start to look like a U-shaped valley. 

In general, different plane wave components of the source will couple to each other through 
Eq. (3). To simulate the gain-compensated lens under more real- world conditions, it was 
therefore tested with several waves traversing the lens simultanously. The transmission of 
one wave as a function of kx-, in the presence of another wave —kx-, is shown in Fig. 2. The 
amplitudes of both waves were set to 1/2 to keep the total field at the source equal to the case 
with a single wave. Moreover, from a number of simulations with several waves, a useful rule 
of thumb was discovered: As a worst-case estimate, one can judge whether the lens operates 
as required by assuming that the mode with largest k^ has amplitude equal to the sum of 
the amplitudes at the source. More preciesly, suppose a single mode kx with amplitude 1 
experiences a transmission greater than 1/2. For any superposition of modes with transversal 
wavenumbers less than kx and sum of amplitudes equal to unity, each mode will experience 
a transmission greater than 1/2. 

For conventional lenses, the Rayleigh criterion is usually applied to quantify the distance 
between two point sources (or in the one-dimensional case, line sources) in order to resolve 
their images. Since our lens is nonlinear, the image of two line sources cannot be determined as 
a superposition of the fields associated with the two sources separately, or as a superposition 
of the fields associated with their Fourier components. Therefore, as in previous literature 
on perfect lenses, we have chosen to consider a single Fourier component source, and defined 
the spatial resolution as 211 /kx-, where kx is the half maximum wavenumber (Fig. 3). Note 
that the behavior of the lens for more complex sources can be determined from the rule of 
thumb described above. 

Fig. 6 shows the absolute value of the transmitted magnetic field at the image plane, 
\H{x,2d)\, resulting from a source consisting of two slits. The image of the slits are clearly 
better resolved with an increased saturation constant Eg. 



4. Conclusion 

We have developed a method for calculating the transmission, reflection, and detailed field 
profile of a gain-compensated perfect lens, taking into account gain saturation. The gain 
compensation clearly improves the resolution limit of perfect lenses. However, due to gain 
saturation, a number of nonideal effects arise, included limited resolution and reflections. 
The nonideal effects depend heavily on the saturation constant and/or the fleld strength of 
the source. 

If there are different waves traversing the lens at the same time, they will interact through 
the material. Waves with a spatial frequency close to the resolution limit will have the 
greatest impact. As a rule of thumb, it is enough to know the sum of amplitudes of the 
waves at the source, and then assume the mode with the largest spatial frequency has this 
amplitude. If this single wave is transmitted, in the sense of a transmission larger than 1/2, 
then so will any superposition of waves with less spatial frequencies and the same sum of 
amplitudes. 

The calculations in this work was performed for TM polarization and a one-dimensional 
source. For a two-dimensional source with both polarizations, both dielectric and magnetic 
losses should be compensated, that is, Im e and Im fi must be reduced. Although the theory in 
this paper can trivially be extended to this situation, there may be serious practical problems 
associated with the fabrication of such active media for optical frequencies. 

For a noncompensated lens, the maximum spatial frequency resolved by the lens is approx- 
imately — (l/d) ln(|l + e|/2) [11,12]. Thus, for a flxed d, an exponential decrease in the losses 
is necessary to increase the resolution linearly. From our numerical results, a similar relation 
is approximately valid for the saturation constant of a gain compensated medium; to achieve 
a linear improvement in the resolution, the saturation constant must increase exponentially. 
This clearly shows the difficulties of achieving very high resolution. 
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List of Figure Captions 

Fig. 1. (Color online) Perfect lens in vacuum. The parameter d is the thickness of the lens, 
a and b are the distances from the source to the lens, and from the lens to the image plane, 
respectively. The parameters are governed by the equation d = a + b. The numbers 1 through 
N indicate the different slices. The lens is considered to be infinite in the xy-plane. 

Fig. 2. (Color online) The absolute value of the transmission coefficient when uq/c = 1 
(normalized), oj^d/c = 27r/10, a = b = d/2, ImXp(a;o) = 0.05, and A^ = 20: (a) Non- 
compensated lens; (6) £s = 10, Ax = 0; (c) £s = 4, Ax = 0; (d) £^ = 10, Ax = 0.015; (e) 
Ss = 10, Ax = 0, two waves, k^ and —kr^, both having amplitude 1/2. 

Fig. 3. (Color online) The resolution of the lens as a function of the saturation constant. 
The resolution is defined as the kx-value where the transmission equals 1/2. Parameters: 
uo/c = 1, uod/c = 27r/10, ImXp(wo) = 0.05, Ax = 0, a = 6 = d/2, and A^ = 20. 

Fig. 4. (Color online) The absolute value of the refiection coefficient at the source plane 
after convergence, for the same cases as those in Fig. 2. 

Fig. 5. (Color online) The distribution of the two components of the evanescent field in 
the lens, for one wave with k^ = 5.1408. The distance is normalized with respect to lens 
thickness, d. Parameters: uq/c = 1, u^d/c = 27r/10, ImXp(i^o) = 0.05, Ax = 0, Ss = 10, 
a = b = d/2, and A^ = 20. The solid line shows the absolute value of the nonzero component 
of h"*" , and the dotted line shows the absolute value of the nonzero component of h~ . 

Fig. 6. (Color online) The absolute value of the transmitted magnetic field at the im- 
age plane, when the source consist of two slits. Parameters: ojq/c = 1, Uod/c = 27r/10, 
ImXp(i^o) = 0.05, Ax = 0, a = b = d/2, and A^ = 20: (a) The incident magnetic field at the 
source, (b) Sg = 0.1, and (c) Sg = 50. 



13 




Fig. 1: (Color online) Perfect lens in vacuum. The parameter d is the thickness of the lens, 
a and h are the distances from the source to the lens, and from the lens to the image 
plane, respectively. The parameters are governed by the equation d = a + h. The numbers 1 
through N indicate the different slices. The lens is considered to be infinite in the x^z-plane. 
ASHKSFl.eps. 
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Fig. 2: (Color online) The absolute value of the transmission coefficient when uq/c = 1 
(normalized), u^d/c = 27r/10, a = b = d/2, ImXp(co'o) = 0.05, and A^ = 20: (a) Non- 
compensated lens; (b) S^ = 10, Ax = 0; (c) S, = 4, Ax = 0; (d) S^ = 10, Ax = 0.015; (e) 
Ss = 10, Ax = 0, two waves, k^ and —kx, both having amplitude 1/2. ASHKSF2.eps. 
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Fig. 3: (Color online) The resolution of the lens as a function of the saturation constant. The 
resolution is defined as the kx-valne where the transmission equals 1/2. Parameters: uq/c = 1, 
uod/c = 27r/10, Imxp(wo) = 0.05, Ax = 0, a = b = d/2, and A^ = 20. ASHKSF3.eps. 
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Fig. 4: (Color online) The absolute value of the reflection coefficient at the source plane after 
convergence, for the same cases as those in Fig. 2. ASHKSF4.eps. 
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Fig. 5: (Color online) The distribution of the two components of the evanescent field in the 
lens, for one wave with kx = 5.1408. The distance is normalized with respect to lens thickness, 
d. Parameters: uq/c = 1, uj^d/c = 27r/10, ImXp(cJo) = 0.05, Ax = 0, £s = 10, a = 6 = d/2, 
and A^ = 20. The solid line shows the absolute value of the nonzero component of h+, and 
the dotted line shows the absolute value of the nonzero component of h~. ASHKSF5.eps. 
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Fig. 6: (Color online) The absolute value of the transmitted magnetic field at the image plane, 
when the source consist of two slits. Parameters: Uq/c = 1, Uod/c = 27r/10, lm.Xp{uJo) = 0.05, 
Ax = 0, a = b = d/2, and A^ = 20: (a) The incident magnetic field at the source, (b) Ss = 0.1, 
and (c) Ss = 50. ASHKSF6.eps. 
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